




%%%%THE FOLLOWING WILL BE USED TO COMPUTE DESCRIPTIVE STATISTICS FOR WEBSITE ACTIVITY:
temp  = g2(~isnan(g2.pre_rawscore) & g2.group2==1 & (g2.contract1==1|g2.contract2==1|g2.contract3==1),:) ;
temp.time_study = (temp.time_study_pre + temp.time_study_post)/2 ;
%%%%FOR SOME REASON, THERE ARE A COUPLE CONTRACT 3 PEOPLE (e.g., ID 2890) WITH 2 PASSED QUIZZES BUT LOGGED TIME OF ZERO.
%%%%WE RE-DEFINE SUCH PEOPLE AS HAVING NO PASSED QUIZZES, RATHER THAN IMPUTING TIME AND DEVICE CONNECTIONS:
yflag = find(temp.total_pass>=2 & temp.timing_pagesubmit==0) ;
temp.total_pass(yflag) = 0 ;
temp.total_attempt(yflag) = 0 ;


DataSize = size(temp)
disp('**********************************************************************')
disp('**** RANK CORRELATIONS BETWEEN SUCCESS RATIO AND MINSPENT: ***********')
disp('**********************************************************************')
[corr_succratlogminspent,p_succratlogminspent] = corr(temp.succrat(temp.total_pass>=2),log(temp.minspent(temp.total_pass>=2)),'type','Spearman','rows','complete')  %#ok
[corr_succrat2minspent,p_succrat2minspent]     = corr(temp.total_pass(temp.pre_rawscore>-inf & temp.total_pass>=2 & temp.group2==1)./...
                                   (temp.minspent(temp.pre_rawscore>-inf & temp.total_pass>=2 & temp.group2==1)),(temp.minspent(temp.total_pass>=2)),'type','Spearman','rows','complete')  %#ok
[corr_succrat2totpass,p_succrat2totpass]       = corr(temp.total_pass(temp.pre_rawscore>-inf & temp.total_pass>=2 & temp.group2==1)./...
                                   (temp.minspent(temp.pre_rawscore>-inf & temp.total_pass>=2 & temp.group2==1)),(temp.total_pass(temp.total_pass>=2)),'type','Spearman','rows','complete')  %#ok


disp('****************************************************************')
disp('********************** DESC STATS: TABLE 1  ********************')
disp('****************************************************************')
stats_active   = [mean(temp.total_pass>=2) median(temp.total_pass>=2) std(temp.total_pass>=2) length(temp.total_pass>=2)...
                  mean(temp.total_pass(temp.contract1==1)>=2)  mean(temp.total_pass(temp.contract2==1)>=2) mean(temp.total_pass(temp.contract3==1)>=2)]  %#ok
stats_onepass  = [mean(temp.total_pass==1) median(temp.total_pass==1) std(temp.total_pass==1) length(temp.total_pass==1)...
                  mean(temp.total_pass(temp.contract1==1)==1)  mean(temp.total_pass(temp.contract2==1)==1) mean(temp.total_pass(temp.contract3==1)==1)]  %#ok
stats_inactive = [mean(temp.total_pass==0) median(temp.total_pass==0) std(temp.total_pass==0) length(temp.total_pass==0)...
                  mean(temp.total_pass(temp.contract1==1)==0)  mean(temp.total_pass(temp.contract2==1)==0) mean(temp.total_pass(temp.contract3==1)==0)]  %#ok
stats_passed   = [mean(temp.total_pass) median(temp.total_pass) std(temp.total_pass) length(temp.total_pass)...
                  mean(temp.total_pass(temp.contract1==1))  mean(temp.total_pass(temp.contract2==1)) mean(temp.total_pass(temp.contract3==1))]  %#ok
disp('%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%')
disp('%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%              
%%%%THE FOLLOWING BLOCK OF CODE IS TO COMPUTE TWO-SAMPLE KS-TESTS FOR DISTRIBUTIONAL EQUALITY:
%%%%Note: We follow the bootstrap procedure of Barrett and Donald (2003, Econometrica, vol. 71(1))
N1 = sum(temp.contract1==1) ;
N2 = sum(temp.contract2==1) ;
N3 = sum(temp.contract3==1) ;
passeddom = (0:1:80)' ;
Fpassed1  = nan(size(passeddom)) ;
Fpassed2  = nan(size(passeddom)) ;
Fpassed3  = nan(size(passeddom)) ;
for ii=1:length(passeddom)
    Fpassed1(ii) = sum( temp.total_pass(temp.contract1==1) <= passeddom(ii) )/N1 ;
    Fpassed2(ii) = sum( temp.total_pass(temp.contract2==1) <= passeddom(ii) )/N2 ;
    Fpassed3(ii) = sum( temp.total_pass(temp.contract3==1) <= passeddom(ii) )/N3 ;
end
S = 10000 ;  %this is the number of bootstrap samples to compute...
%%%% END OF KS TEST SETUP CODE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PassedContract12_diff = mean(temp.total_pass(temp.contract2==1)) - mean(temp.total_pass(temp.contract1==1))  %#ok
PassedContract12_pctdiff = 100*(mean(temp.total_pass(temp.contract2==1)) - mean(temp.total_pass(temp.contract1==1)))/mean(temp.total_pass(temp.contract1==1))  %#ok
PassedContract12_pooled_SE = sqrt( var(temp.total_pass(temp.contract2==1))/(N2) + ...
                          var(temp.total_pass(temp.contract1==1))/(N1) )   %#ok
PassedContract12TStat = ( mean(temp.total_pass(temp.contract2==1)) - mean(temp.total_pass(temp.contract1==1)) )/PassedContract12_pooled_SE  %#ok
PassedContract12PValue = 2*normcdf(-abs(PassedContract12TStat))  %#ok
tic
disp('%%%%BARRETT-DONALD TEST FOR FIRST-ORDER DOMINANCE OF QUIZ OUTPUT IN CONTRACT 2 VERSUS CONTRACT 1:')
PassedContract12KSBDStat = max(Fpassed1 - Fpassed2)*sqrt( (N1*N2)/(N1+N2) )  %#ok
tempsamp = temp.total_pass(temp.contract1==1 | temp.contract2==1) ; 
Sbar12 = nan(S,1) ; 
for s=1:S
    temp1 = datasample(tempsamp,N1) ; 
    temp2 = datasample(tempsamp,N2) ; 
    Ftemp1 = nan(size(passeddom)) ; 
    Ftemp2 = nan(size(passeddom)) ; 
    for ii=1:length(passeddom) 
        Ftemp1(ii) = sum( temp1 <= passeddom(ii) )/N1 ; 
        Ftemp2(ii) = sum( temp2 <= passeddom(ii) )/N2 ; 
    end 
    Sbar12(s) = max(Ftemp1 - Ftemp2)*sqrt( (N1*N2)/(N1+N2) ) ;
end 
PassedContract12KSPval = mean( Sbar12 > PassedContract12KSBDStat )  %#ok 
Nbootstraps12 = S %#ok
toc

PassedContract13_diff = mean(temp.total_pass(temp.contract3==1)) - mean(temp.total_pass(temp.contract1==1))  %#ok
PassedContract13_pctdiff = 100*(mean(temp.total_pass(temp.contract3==1)) - mean(temp.total_pass(temp.contract1==1)))/mean(temp.total_pass(temp.contract1==1))  %#ok
PassedContract13_pooled_SE = sqrt( var(temp.total_pass(temp.contract3==1))/(N3) + ...
                          var(temp.total_pass(temp.contract1==1))/(N1) )   %#ok
PassedContract13TStat = ( mean(temp.total_pass(temp.contract3==1)) - mean(temp.total_pass(temp.contract1==1)) )/PassedContract13_pooled_SE  %#ok
PassedContract13PValue = 2*normcdf(-abs(PassedContract13TStat))  %#ok
tic
disp('%%%%BARRETT-DONALD TEST FOR FIRST-ORDER DOMINANCE OF QUIZ OUTPUT IN CONTRACT 3 VERSUS CONTRACT 1:')
PassedContract13KSBDStat = max(Fpassed1 - Fpassed3)*sqrt( (N1*N3)/(N1+N3) )  %#ok
tempsamp = temp.total_pass(temp.contract1==1 | temp.contract3==1) ; 
Sbar13 = nan(S,1) ; 
for s=1:S
    temp1 = datasample(tempsamp,N1) ; 
    temp3 = datasample(tempsamp,N3) ; 
    Ftemp1 = nan(size(passeddom)) ; 
    Ftemp3 = nan(size(passeddom)) ; 
    for ii=1:length(passeddom) 
        Ftemp1(ii) = sum( temp1 <= passeddom(ii) )/N1 ; 
        Ftemp3(ii) = sum( temp3 <= passeddom(ii) )/N3 ; 
    end 
    Sbar13(s) = max(Ftemp1 - Ftemp3)*sqrt( (N1*N3)/(N1+N3) ) ;
end 
PassedContract13KSPval = mean( Sbar13 > PassedContract13KSBDStat )  %#ok 
Nbootstraps13 = S %#ok
toc

PassedContract23_diff = mean(temp.total_pass(temp.contract3==1)) - mean(temp.total_pass(temp.contract2==1))  %#ok
PassedContract23_pctdiff = 100*(mean(temp.total_pass(temp.contract3==1)) - mean(temp.total_pass(temp.contract2==1)))/mean(temp.total_pass(temp.contract2==1))  %#ok
PassedContract23_pooled_SE = sqrt( var(temp.total_pass(temp.contract3==1))/(N3) + ...
                          var(temp.total_pass(temp.contract2==1))/(N2) )   %#ok
PassedContract23TStat = ( mean(temp.total_pass(temp.contract3==1)) - mean(temp.total_pass(temp.contract2==1)) )/PassedContract23_pooled_SE  %#ok
PassedContract23PValue = 2*normcdf(-abs(PassedContract23TStat))  %#ok
tic
disp('%%%%BARRETT-DONALD TEST FOR FIRST-ORDER DOMINANCE OF QUIZ OUTPUT IN CONTRACT 3 VERSUS CONTRACT 2:')
PassedContract23KSBDStat = max(Fpassed2 - Fpassed3)*sqrt( (N2*N3)/(N2+N3) )  %#ok
tempsamp = temp.total_pass(temp.contract2==1 | temp.contract3==1) ; 
Sbar23 = nan(S,1) ; 
for s=1:S
    temp2 = datasample(tempsamp,N2) ; 
    temp3 = datasample(tempsamp,N3) ; 
    Ftemp2 = nan(size(passeddom)) ; 
    Ftemp3 = nan(size(passeddom)) ; 
    for ii=1:length(passeddom) 
        Ftemp2(ii) = sum( temp2 <= passeddom(ii) )/N2 ; 
        Ftemp3(ii) = sum( temp3 <= passeddom(ii) )/N3 ; 
    end 
    Sbar23(s) = max(Ftemp2 - Ftemp3)*sqrt( (N2*N3)/(N2+N3) ) ;
end 
PassedContract23KSPval = mean( Sbar23 > PassedContract23KSBDStat )  %#ok 
Nbootstraps23 = S %#ok
toc
disp('%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%')
disp('%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%')
disp('\n')
disp('\n')
disp('\n')
disp('\n')


disp('%%%%%%%%%% NOW RE-DO ALL THAT KS TEST STUFF USING THE CONTRACT-3-EQUIVALENT SET OF UPTAKERS...')
disp('%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%')
disp('%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%')
refmass = stats_active(end)  %#ok
samptrunc2 = round(N2*refmass)  %#ok
samptrunc1 = round(N1*refmass)  %#ok
tempC3E_3 = temp.total_pass(temp.contract3==1 & temp.total_pass>=2) ;
tempC3E_2 = maxk(temp.total_pass(temp.contract2==1),samptrunc2) ;
tempC3E_1 = maxk(temp.total_pass(temp.contract1==1),samptrunc1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%              
%%%%THE FOLLOWING BLOCK OF CODE IS TO COMPUTE TWO-SAMPLE KS-TESTS FOR DISTRIBUTIONAL EQUALITY:
%%%%Note: We follow the bootstrap procedure of Barrett and Donald (2003, Econometrica, vol. 71(1))
N1 = length(tempC3E_1) ;
N2 = length(tempC3E_2) ;
N3 = length(tempC3E_3) ;
passeddom = (0:1:80)' ;
Fpassed1  = nan(size(passeddom)) ;
Fpassed2  = nan(size(passeddom)) ;
Fpassed3  = nan(size(passeddom)) ;
for ii=1:length(passeddom)
    Fpassed1(ii) = sum( tempC3E_1 <= passeddom(ii) )/N1 ;
    Fpassed2(ii) = sum( tempC3E_2 <= passeddom(ii) )/N2 ;
    Fpassed3(ii) = sum( tempC3E_3 <= passeddom(ii) )/N3 ;
end
S = 100000 ;  %this is the number of bootstrap samples to compute...

%%%% END OF KS TEST SETUP CODE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PassedContract12_diffC3E = mean(tempC3E_2) - mean(tempC3E_1)  %#ok
PassedContract12_pctdiffC3E = 100*(mean(tempC3E_2) - mean(tempC3E_1))/mean(tempC3E_1)  %#ok
PassedContract12_pooled_SEC3E = sqrt( var(tempC3E_2)/(N2) + ...
                          var(tempC3E_1)/(N1) )   %#ok
PassedContract12TStatC3E = ( mean(tempC3E_2) - mean(tempC3E_1) )/PassedContract12_pooled_SEC3E  %#ok
PassedContract12PValueC3E = 2*normcdf(-abs(PassedContract12TStatC3E))  %#ok
tic
disp('%%%%BARRETT-DONALD TEST FOR FIRST-ORDER DOMINANCE OF QUIZ OUTPUT IN CONTRACT 2 VERSUS CONTRACT 1:')
PassedContract12KSBDStatC3E = max(Fpassed1 - Fpassed2)*sqrt( (N1*N2)/(N1+N2) )  %#ok
tempsamp = [tempC3E_1; tempC3E_2] ; 
Sbar12C3E = nan(S,1) ; 
for s=1:S
    temp1 = datasample(tempsamp,N1) ; 
    temp2 = datasample(tempsamp,N2) ; 
    Ftemp1 = nan(size(passeddom)) ; 
    Ftemp2 = nan(size(passeddom)) ; 
    for ii=1:length(passeddom) 
        Ftemp1(ii) = sum( temp1 <= passeddom(ii) )/N1 ; 
        Ftemp2(ii) = sum( temp2 <= passeddom(ii) )/N2 ; 
    end 
    Sbar12C3E(s) = max(Ftemp1 - Ftemp2)*sqrt( (N1*N2)/(N1+N2) ) ;
end 
PassedContract12KSPvalC3E = mean( Sbar12C3E > PassedContract12KSBDStatC3E )  %#ok 
Nbootstraps12C3E = S %#ok 
toc

PassedContract13_diffC3E = mean(tempC3E_3) - mean(tempC3E_1)  %#ok 
PassedContract13_pctdiffC3E = 100*(mean(tempC3E_3) - mean(tempC3E_1))/mean(tempC3E_1)  %#ok 
PassedContract13_pooled_SEC3E = sqrt( var(tempC3E_3)/(N3) + ...
                          var(tempC3E_1)/(N1) )   %#ok 
PassedContract13TStatC3E = ( mean(tempC3E_3) - mean(tempC3E_1) )/PassedContract13_pooled_SEC3E  %#ok 
PassedContract13PValueC3E = 2*normcdf(-abs(PassedContract13TStatC3E))  %#ok 
tic
disp('%%%%BARRETT-DONALD TEST FOR FIRST-ORDER DOMINANCE OF QUIZ OUTPUT IN CONTRACT 3 VERSUS CONTRACT 1:')
PassedContract13KSBDStatC3E = max(Fpassed1 - Fpassed3)*sqrt( (N1*N3)/(N1+N3) )  %#ok 
tempsamp = [tempC3E_1; tempC3E_3] ; 
Sbar13C3E = nan(S,1) ; 
for s=1:S
    temp1 = datasample(tempsamp,N1) ; 
    temp3 = datasample(tempsamp,N3) ; 
    Ftemp1 = nan(size(passeddom)) ; 
    Ftemp3 = nan(size(passeddom)) ; 
    for ii=1:length(passeddom) 
        Ftemp1(ii) = sum( temp1 <= passeddom(ii) )/N1 ; 
        Ftemp3(ii) = sum( temp3 <= passeddom(ii) )/N3 ; 
    end 
    Sbar13C3E(s) = max(Ftemp1 - Ftemp3)*sqrt( (N1*N3)/(N1+N3) ) ;
end 
PassedContract13KSPvalC3E = mean( Sbar13C3E > PassedContract13KSBDStatC3E )  %#ok 
Nbootstraps13C3E = S %#ok
toc

PassedContract23_diffC3E = mean(tempC3E_3) - mean(tempC3E_2)  %#ok
PassedContract23_pctdiffC3E = 100*(mean(tempC3E_3) - mean(tempC3E_2))/mean(tempC3E_2)  %#ok
PassedContract23_pooled_SEC3E = sqrt( var(tempC3E_3)/(N3) + ...
                          var(tempC3E_2)/(N2) )   %#ok
PassedContract23TStatC3E = ( mean(tempC3E_3) - mean(tempC3E_2) )/PassedContract23_pooled_SEC3E  %#ok
PassedContract23PValueC3E = 2*normcdf(-abs(PassedContract23TStatC3E))  %#ok
tic
disp('%%%%BARRETT-DONALD TEST FOR FIRST-ORDER DOMINANCE OF QUIZ OUTPUT IN CONTRACT 3 VERSUS CONTRACT 2:')
PassedContract23KSBDStatC3E = max(Fpassed2 - Fpassed3)*sqrt( (N2*N3)/(N2+N3) )  %#ok
tempsamp = [tempC3E_3; tempC3E_2] ; 
Sbar23C3E = nan(S,1) ; 
for s=1:S
    temp2 = datasample(tempsamp,N2) ; 
    temp3 = datasample(tempsamp,N3) ; 
    Ftemp2 = nan(size(passeddom)) ; 
    Ftemp3 = nan(size(passeddom)) ; 
    for ii=1:length(passeddom) 
        Ftemp2(ii) = sum( temp2 <= passeddom(ii) )/N2 ; 
        Ftemp3(ii) = sum( temp3 <= passeddom(ii) )/N3 ; 
    end 
    Sbar23C3E(s) = max(Ftemp2 - Ftemp3)*sqrt( (N2*N3)/(N2+N3) ) ;
end 
PassedContract23KSPvalC3E = mean( Sbar23C3E > PassedContract23KSBDStatC3E )  %#ok 
Nbootstraps23C3E = S %#ok
toc
disp('**************************************************************************') 
disp('***** FINALLY, RE-DO ACTIVE ONLY STATS FOR CONTRACT 3 EQUIVALENT SAMPLE...') 
disp('**************************************************************************') 
stats_passedACTC3E   = [mean([tempC3E_1; tempC3E_2; tempC3E_3]) median([tempC3E_1; tempC3E_2; tempC3E_3])...
                        std([tempC3E_1; tempC3E_2; tempC3E_3]) length([tempC3E_1; tempC3E_2; tempC3E_3])...
                        mean(tempC3E_1)  mean(tempC3E_2) mean(tempC3E_3)]  %#ok  



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%               FIGURE 1:                   %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
 box on; set(gca,'xgrid','on'); set(gca,'ygrid','on')
hold on
[F_PASS_1,X_PASS_1] = ecdf(tempC3E_1) ; 
[F_PASS_2,X_PASS_2] = ecdf(tempC3E_2) ; 
[F_PASS_3,X_PASS_3] = ecdf(tempC3E_3) ; 
stairs(X_PASS_1,F_PASS_1,'-.k','LineWidth',3)
stairs(X_PASS_2,F_PASS_2,'b','LineWidth',3)
stairs(X_PASS_3,F_PASS_3,'--r','LineWidth',3)
ylim([0,1])
xlabel('TOTAL COMPLETED LEARNING TASKS','FontSize',20,'Fontweight','bold')
ylabel('CONDITIONAL CDFs','FontSize',20,'Fontweight','bold')
legend('CONTRACT GROUP 1','CONTRACT GROUP 2','CONTRACT GROUP 3','FontSize',20,'FontWeight','bold','Location','SouthEast')

%%
disp('%%%%%  FINALLY, FOR ILLUSTRATIVE PURPOSES, HERE ARE THE INTEGRAL DIFFERENCES FOR THE QUANTILE FUNCTIONS:  %%%%%')
commondom12 = unique([F_PASS_1;F_PASS_2]); 
Qdiff12 = interp1(F_PASS_2,X_PASS_2,commondom12,'previous') - interp1(F_PASS_1,X_PASS_1,commondom12,'previous') ;
QuantileDiffIntACTC3E12 = sum( (commondom12(2:end)-commondom12(1:end-1)).*Qdiff12(1:end-1) ) %#ok  %NOTE: REIMANN INTEGRATION IS FINE HERE BECAUSE WE ARE WORKING WITH STEP FUNCTIONS...
commondom13 = unique([F_PASS_1;F_PASS_3]); 
Qdiff13 = interp1(F_PASS_3,X_PASS_3,commondom13,'previous') - interp1(F_PASS_1,X_PASS_1,commondom13,'previous') ;
QuantileDiffIntACTC3E13 = sum( (commondom13(2:end)-commondom13(1:end-1)).*Qdiff13(1:end-1) ) %#ok  %NOTE: REIMANN INTEGRATION IS FINE HERE BECAUSE WE ARE WORKING WITH STEP FUNCTIONS...
commondom23 = unique([F_PASS_2;F_PASS_3]); 
Qdiff23 = interp1(F_PASS_3,X_PASS_3,commondom23,'previous') - interp1(F_PASS_2,X_PASS_2,commondom23,'previous') ;
QuantileDiffIntACTC3E23 = sum( (commondom23(2:end)-commondom23(1:end-1)).*Qdiff23(1:end-1) ) %#ok  %NOTE: REIMANN INTEGRATION IS FINE HERE BECAUSE WE ARE WORKING WITH STEP FUNCTIONS...
SECOND_QDiffIntACTC3E_12versus23 = QuantileDiffIntACTC3E23 - QuantileDiffIntACTC3E12 %#ok
QDiffPctChangeACTC3E = 100*SECOND_QDiffIntACTC3E_12versus23/QuantileDiffIntACTC3E12 %#ok
disp('%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%')
disp('%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%BARRETT-DONALDEcma2003KSTest%%%%')




              
stats_solved   = [mean(6*temp.total_pass) median(6*temp.total_pass) std(6*temp.total_pass) length(6*temp.total_pass)...
                  mean(6*temp.total_pass(temp.contract1==1))  mean(6*temp.total_pass(temp.contract2==1)) mean(6*temp.total_pass(temp.contract3==1))]  %#ok
              
stats_minspent = [mean(temp.minspent) median(temp.minspent) std(temp.minspent) length(temp.minspent)...
                  mean(temp.minspent(temp.contract1==1))  mean(temp.minspent(temp.contract2==1)) mean(temp.minspent(temp.contract3==1))]  %#ok  
              
temp.time_study = (temp.time_study_pre+temp.time_study_post)/2 ;              
stats_hwkALL   = [mean(temp.time_study) median(temp.time_study) std(temp.time_study) length(temp.time_study)...
                  mean(temp.time_study(temp.contract1==1))  mean(temp.time_study(temp.contract2==1)) mean(temp.time_study(temp.contract3==1))]  %#ok
              
stats_hwkALLpre= [mean(temp.time_study_pre) median(temp.time_study_pre) std(temp.time_study_pre) length(temp.time_study_pre)...
                  mean(temp.time_study_pre(temp.contract1==1))  mean(temp.time_study_pre(temp.contract2==1)) mean(temp.time_study_pre(temp.contract3==1))]  %#ok
              
stats_hwkALLpost= [mean(temp.time_study_post) median(temp.time_study_post) std(temp.time_study_post) length(temp.time_study_post)...
                  mean(temp.time_study_post(temp.contract1==1))  mean(temp.time_study_post(temp.contract2==1)) mean(temp.time_study_post(temp.contract3==1))]  %#ok
              
hwkprepostALL_diff = ( mean(temp.time_study_pre) - mean(temp.time_study_post) )  %#ok

hwkprepostALL_pctdiff = 100*( mean(temp.time_study_pre) - mean(temp.time_study_post) )/mean(temp.time_study_pre)  %#ok

hwkprepostALL_pooled_SE = sqrt( var(temp.time_study_pre)/(length(temp.time_study_pre)) + ...
                          var(temp.time_study_post)/(length(temp.time_study_post)) )   %#ok
                      
PrePostHwkALLTStat = ( mean(temp.time_study_pre) - mean(temp.time_study_post) )/hwkprepostALL_pooled_SE  %#ok

PrePostHwkALLPValue = 2*normcdf(-abs(PrePostHwkALLTStat))  %#ok

stats_hwkACT   = [mean(temp.time_study(temp.total_pass>=2)) median(temp.time_study(temp.total_pass>=2)) std(temp.time_study(temp.total_pass>=2)) length(temp.time_study(temp.total_pass>=2))...
                  mean(temp.time_study(temp.total_pass>=2 & temp.contract1==1))  mean(temp.time_study(temp.total_pass>=2 & temp.contract2==1)) mean(temp.time_study(temp.total_pass>=2 & temp.contract3==1))]  %#ok

stats_hwkpreACT= [mean(temp.time_study_pre(temp.total_pass>=2)) median(temp.time_study_pre(temp.total_pass>=2)) std(temp.time_study_pre(temp.total_pass>=2)) length(temp.time_study_pre(temp.total_pass>=2))...
                  mean(temp.time_study_pre(temp.total_pass>=2 & temp.contract1==1))  mean(temp.time_study_pre(temp.total_pass>=2 & temp.contract2==1)) mean(temp.time_study_pre(temp.total_pass>=2 & temp.contract3==1))]  %#ok
              
stats_hwkpostACT= [mean(temp.time_study_post(temp.total_pass>=2)) median(temp.time_study_post(temp.total_pass>=2)) std(temp.time_study_post(temp.total_pass>=2)) length(temp.time_study_post(temp.total_pass>=2))...
                  mean(temp.time_study_post(temp.total_pass>=2 & temp.contract1==1))  mean(temp.time_study_post(temp.total_pass>=2 & temp.contract2==1)) mean(temp.time_study_post(temp.total_pass>=2 & temp.contract3==1))]  %#ok
              
hwkprepostACT_diff = ( mean(temp.time_study_pre(temp.total_pass>=2)) - mean(temp.time_study_post(temp.total_pass>=2)) )  %#ok

hwkprepostACT_pctdiff = 100*( mean(temp.time_study_pre(temp.total_pass>=2)) - mean(temp.time_study_post(temp.total_pass>=2)) )/mean(temp.time_study_pre(temp.total_pass>=2))  %#ok

hwkprepostACT_pooled_SE = sqrt( var(temp.time_study_pre(temp.total_pass>=2))/(length(temp.time_study_pre(temp.total_pass>=2))) + ...
                          var(temp.time_study_post(temp.total_pass>=2))/(length(temp.time_study_post(temp.total_pass>=2))) )   %#ok
                      
PrePostHwkACTTStat = ( mean(temp.time_study_pre(temp.total_pass>=2)) - mean(temp.time_study_post(temp.total_pass>=2)) )/hwkprepostACT_pooled_SE    %#ok

PrePostHwkACTPValue = 2*normcdf(-abs(PrePostHwkACTTStat))   %#ok

stats_hwk1Pass = [mean(temp.time_study(temp.total_pass==1)) median(temp.time_study(temp.total_pass==1)) std(temp.time_study(temp.total_pass==1)) length(temp.time_study(temp.total_pass==1))...
                  mean(temp.time_study(temp.total_pass==1 & temp.contract1==1))  mean(temp.time_study(temp.total_pass==1 & temp.contract2==1)) mean(temp.time_study(temp.total_pass==1 & temp.contract3==1))]  %#ok
              
stats_hwkINkACT= [mean(temp.time_study(temp.total_pass<2)) median(temp.time_study(temp.total_pass<2)) std(temp.time_study(temp.total_pass<2)) length(temp.time_study(temp.total_pass<2))...
                  mean(temp.time_study(temp.total_pass<2 & temp.contract1==1))  mean(temp.time_study(temp.total_pass<2 & temp.contract2==1)) mean(temp.time_study(temp.total_pass<2 & temp.contract3==1))]   %#ok          
              
hwkActInact_diff       = mean(temp.time_study(temp.total_pass>=2)) - mean(temp.time_study(temp.total_pass<2))  %#ok

hwkActInact_pctdiff   = 100*( mean(temp.time_study(temp.total_pass>=2)) - mean(temp.time_study(temp.total_pass<2)) )/mean(temp.time_study(temp.total_pass>=2))  %#ok

hwkActInact_pooled_SE = sqrt( var(temp.time_study(temp.total_pass>=2))/(length(temp.time_study(temp.total_pass>=2))) + ...
                          var(temp.time_study(temp.total_pass<2))/(length(temp.time_study(temp.total_pass<2))-1) )   %#ok
                      
HwkActInactTStat = ( mean(temp.time_study(temp.total_pass>=2)) - mean(temp.time_study_post(temp.total_pass<2)) )/hwkActInact_pooled_SE    %#ok

HwkActInactPValue = 2*normcdf(-abs(HwkActInactTStat))   %#ok




disp('***************************************************************************************')
disp('*************** CORRELATIONS BETWEEN HOMEWORK TIME AND MATH ACTIVITY:  ****************')
disp('***************************************************************************************')
[CorrHwkTotPassALL,PvalHwkTotPassALL]                 = corr(temp.time_study,temp.total_pass,'Type','Pearson','Rows','complete')  %#ok

[CorrHwkTotPassALLSpearman,PvalHwkTotPassALLSpearman] = corr(temp.time_study,temp.total_pass,'Type','Spearman','Rows','complete')  %#ok

[CorrHwkMinspentALL,PvalHwkMinspentALL]                 = corr(temp.time_study,temp.minspent,'Type','Pearson','Rows','complete')  %#ok

[CorrHwkMinspentALLSpearman,PvalHwkMinspentALLSpearman] = corr(temp.time_study,temp.minspent,'Type','Spearman','Rows','complete')  %#ok

[CorrFreqMissHwkMinspentALL,PvalFreqMissHwkMinspentALL]                 = corr(temp.freq_miss_hmwk,temp.minspent,'Type','Pearson','Rows','complete')  %#ok

[CorrFreqMissHwkMinspentALLSpearman,PvalFreqMissHwkMinspentALLSpearman] = corr(temp.freq_miss_hmwk,temp.minspent,'Type','Spearman','Rows','complete')  %#ok


[CorrHwkActStatus,PvalHwkActStatus]                 = corr(temp.time_study,temp.total_pass>=2,'Type','Pearson','Rows','complete')  %#ok

disp('***************************************************************************************')
disp('***************************************************************************************')
disp('')
disp('')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%               TABLE OS.2:                 %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
temp.mobilefrac   = temp.mobile./(temp.mobile+temp.computer+temp.tablet) ; 
temp.computerfrac = temp.computer./(temp.mobile+temp.computer+temp.tablet) ; 
temp.tabletfrac   = temp.tablet./(temp.mobile+temp.computer+temp.tablet) ; 
disp('***************************************************************************************')
disp('********** EXPLORATORY ANALYSIS ON THOSE WHO REPORTED NO INTERNET AT HOME:  ***********')
disp('***************************************************************************************')
N_no_internet = sum(temp.no_internet,'omitnan')  %#ok
NO_INTERNET_MASS = N_no_internet/length(temp.student_id)   %#ok
N_participants_no_internet = sum(temp.no_internet==1 & temp.timing_pagesubmit>0)  %#ok
N_participants_no_internet/N_no_internet  %#ok
STATS_no_internetACT = [mean(temp.total_pass(temp.no_internet==1 & temp.timing_pagesubmit>0)) ...
median(temp.total_pass(temp.no_internet==1 & temp.timing_pagesubmit>0)) ...
std(temp.total_pass(temp.no_internet==1 & temp.timing_pagesubmit>0)) ...
length(temp.total_pass(temp.no_internet==1 & temp.timing_pagesubmit>0))]    %#ok
No_internet_computerfrac = mean(temp.computerfrac(temp.no_internet==1 & temp.timing_pagesubmit>0))  %#ok
No_internet_tabletfrac = mean(temp.tabletfrac(temp.no_internet==1 & temp.timing_pagesubmit>0))  %#ok
No_internet_mobilefrac = mean(temp.mobilefrac(temp.no_internet==1 & temp.timing_pagesubmit>0))  %#ok

flag = find(~isnan(temp.no_internet)) ;
Xtemp = [ones(length(temp.no_internet(flag)),1), temp.no_internet(flag), 1*(temp.district_id(flag)==3), 1*(temp.district_id(flag)==4), ...
         temp.mean_inc_cbg(flag), temp.time_study(flag), temp.math_attitude(flag), temp.contract2(flag), temp.contract3(flag)] ;

Ytemp = temp.total_pass(flag) ;
disp('REGRESSION OF total_pass ON NO_INTERNET, DISTRICT DUMMIES, NBHD INCOME, TIME_STUDY, AND MATH_ATTITUDE,')
disp(' AND CONTRACT INCENTIVES. SECOND COEFFICIENT IS NO_INTERNET. WE ALSO COMPUTE STANDARD ERRORS AND P-VALUES.')
[Btotal_pass,Btotal_passINT,~,~,STATStotal_pass] = regress(Ytemp,Xtemp,0.05)  %#ok
SE = (Btotal_passINT(:,2) - Btotal_pass)/abs(norminv(0.025))  %#ok
T  = Btotal_pass./SE  %#ok
Pval = 2*normcdf(-1*abs(T))   %#ok

Xtemp = [ones(length(temp.no_internet(flag)),1), temp.no_internet(flag), 1*(temp.district_id(flag)==3), 1*(temp.district_id(flag)==4), ...
         temp.mean_inc_cbg(flag), temp.mean_inc_cbg(flag).^2, temp.time_study(flag), temp.time_study(flag).^2, temp.math_attitude(flag), temp.math_attitude(flag).^2,...
         temp.contract2(flag), temp.contract3(flag), temp.female(flag), temp.black(flag), temp.hispanic(flag)] ;
Ytemp = temp.total_pass(flag) ;
disp('REGRESSION OF total_pass ON NO_INTERNET, DISTRICT DUMMIES, NBHD INCOME (QUADRATIC),')
disp('TIME_STUDY (QUADRATIC), MATH_ATTITUDE (QUADRATIC), CONTRACT INCENTIVES, GENDER, AND RACE. SECOND')
disp('COEFFICIENT IS NO_INTERNET.  WE ALSO COMPUTE STANDARD ERRORS AND P-VALUES.')
[Btotal_pass,Btotal_passINT,~,~,STATStotal_pass] = regress(Ytemp,Xtemp,0.05)  %#ok
SE = (Btotal_passINT(:,2) - Btotal_pass)/abs(norminv(0.025))  %#ok
T  = Btotal_pass./SE  %#ok
Pval = 2*normcdf(-1*abs(T))    %#ok
disp('***************************************************************************************')
disp('***************************************************************************************')
disp('**************************  END EXPLORATORY ANALYSIS...  ******************************')
disp('***************************************************************************************')
disp('***************************************************************************************')


temp.pay   = zeros(length(temp.total_pass),1);
temp.payPR = zeros(length(temp.total_pass),1); 
for ii=1:length(temp.pay)
    if temp.total_pass(ii)>=2 && temp.contract1(ii)==1
        temp.pay(ii) = 15 + 0.75*temp.total_pass(ii) ;
        temp.payPR(ii) = 0.75*temp.total_pass(ii) ;
    elseif temp.total_pass(ii)>=2 && temp.contract2(ii)==1
        temp.pay(ii) = 10 + 1.00*temp.total_pass(ii) ;
        temp.payPR(ii) = 1.00*temp.total_pass(ii) ;
    elseif temp.total_pass(ii)>=2 && temp.contract3(ii)==1
        temp.pay(ii) = 5 + 1.25*temp.total_pass(ii) ;
        temp.payPR(ii) = 1.25*temp.total_pass(ii) ;
    end
end
stats_payALL      = [mean(temp.pay) median(temp.pay) std(temp.pay) length(temp.pay)...
                  mean(temp.pay(temp.contract1==1))  mean(temp.pay(temp.contract2==1)) mean(temp.pay(temp.contract3==1))]  %#ok

marginal = temp(temp.total_pass==1 & (temp.contract1==1 | temp.contract2==1),:) ; 
disp('************************************************************************') 
disp('***** NOW, DROP NON-ACTIVES AND RE-COMPUTE DESC-STATS FOR ACTIVE ONLY...') 
disp('************************************************************************') 
temp(temp.total_pass<2,:) = [] ; 
stats_passedACT   = [mean(temp.total_pass) median(temp.total_pass) std(temp.total_pass) length(temp.total_pass)...
                     mean(temp.total_pass(temp.contract1==1))  mean(temp.total_pass(temp.contract2==1)) mean(temp.total_pass(temp.contract3==1))]  %#ok
stats_attemptsACT = [mean(temp.total_attempt) median(temp.total_attempt) std(temp.total_attempt) length(temp.total_attempt)...
                     mean(temp.total_attempt(temp.contract1==1))  mean(temp.total_attempt(temp.contract2==1)) mean(temp.total_attempt(temp.contract3==1))]  %#ok
stats_solvedACT   = [mean(6*temp.total_pass) median(6*temp.total_pass) std(6*temp.total_pass) length(6*temp.total_pass)...
                     mean(6*temp.total_pass(temp.contract1==1))  mean(6*temp.total_pass(temp.contract2==1)) mean(6*temp.total_pass(temp.contract3==1))]  %#ok
stats_minspentACT = [mean(temp.minspent) median(temp.minspent) std(temp.minspent) length(temp.minspent)...
                     mean(temp.minspent(temp.contract1==1))  mean(temp.minspent(temp.contract2==1)) mean(temp.minspent(temp.contract3==1))]  %#ok
stats_minpercomp  = [mean((temp.minspent)./temp.total_pass) median((temp.minspent)./temp.total_pass) std((temp.minspent)./temp.total_pass) length((temp.minspent)./temp.total_pass)...
                     mean((temp.minspent(temp.contract1==1))./temp.total_pass(temp.contract1==1))  mean((temp.minspent(temp.contract2==1))./temp.total_pass(temp.contract2==1)) mean((temp.minspent(temp.contract3==1))./temp.total_pass(temp.contract3==1))]  %#ok
stats_payACT      = [mean(temp.pay) median(temp.pay) std(temp.pay) length(temp.pay)...
                     mean(temp.pay(temp.contract1==1))  mean(temp.pay(temp.contract2==1)) mean(temp.pay(temp.contract3==1))]  %#ok
stats_payphrACT   = [mean(temp.pay./(temp.minspent/60)) median(temp.pay./(temp.minspent/60)) std(temp.pay./(temp.minspent/60)) length(temp.pay./(temp.minspent/60))...
                     mean(temp.pay(temp.contract1==1)./(temp.minspent(temp.contract1==1)/60))  mean(temp.pay(temp.contract2==1)./(temp.minspent(temp.contract2==1)/60)) mean(temp.pay(temp.contract3==1)./(temp.minspent(temp.contract3==1)/60))]  %#ok
stats_payPRphrACT = [mean(temp.payPR./(temp.minspent/60)) median(temp.payPR./(temp.minspent/60)) std(temp.payPR./(temp.minspent/60)) length(temp.payPR./(temp.minspent/60))...
                     mean(temp.payPR(temp.contract1==1)./(temp.minspent(temp.contract1==1)/60))  mean(temp.payPR(temp.contract2==1)./(temp.minspent(temp.contract2==1)/60)) mean(temp.payPR(temp.contract3==1)./(temp.minspent(temp.contract3==1)/60))]  %#ok
AllConnections_PCfrac     = sum(temp.computer)/(sum(temp.computer + temp.mobile + temp.tablet))   %#ok
AllConnections_Tabletfrac = sum(temp.tablet)/(sum(temp.computer + temp.mobile + temp.tablet))   %#ok
AllConnections_Mobilefrac = sum(temp.mobile)/(sum(temp.computer + temp.mobile + temp.tablet))   %#ok
stats_PCfrac      = [mean(temp.computerfrac) median(temp.computerfrac) std(temp.computerfrac) length(temp.computerfrac) ...
                     mean(temp.computerfrac(temp.contract1==1)) mean(temp.computerfrac(temp.contract2==1)) mean(temp.computerfrac(temp.contract3==1))]   %#ok
stats_tabletfrac  = [mean(temp.tabletfrac) median(temp.tabletfrac) std(temp.tabletfrac) length(temp.tabletfrac) ...
                     mean(temp.tabletfrac(temp.contract1==1)) mean(temp.tabletfrac(temp.contract2==1)) mean(temp.tabletfrac(temp.contract3==1))]   %#ok
stats_mobilefrac  = [mean(temp.mobilefrac) median(temp.mobilefrac) std(temp.mobilefrac) length(temp.mobilefrac) ...
                     mean(temp.mobilefrac(temp.contract1==1)) mean(temp.mobilefrac(temp.contract2==1)) mean(temp.mobilefrac(temp.contract3==1))]   %#ok

